Fox's Git Mirrors
Commit 33052aae8df1a0bff5422067da6d041c3d038841
Parents : f3d3afb
Author : Ahmet Inan <inan@aicodix.de>
Date : 2024-02-03T10:29:54+01:00
experimenting with list decoding
Changes
2 files changed, 343 insertions(+), 0 deletions(-)
Diff
diff --git a/osd.hh b/osd.hh
index 9548cf2..91b9dcf 100644
--- a/osd.hh
+++ b/osd.hh
@@ -237,5 +237,158 @@ public:
}
};
+template <int N, int K, int O, int L>
+class OrderedStatisticsListDecoder
+{
+ static const int S = sizeof(size_t);
+ static const int W = (N+S-1) & ~(S-1);
+ int8_t G[W*K];
+ int8_t codeword[W], candidate[W*2*L];
+ int8_t softperm[W];
+ int16_t perm[W];
+ int score[2*L], scoreperm[2*L];
+ void row_echelon()
+ {
+ for (int k = 0; k < K; ++k) {
+ // find pivot in this column
+ for (int j = k; j < K; ++j) {
+ if (G[W*j+k]) {
+ for (int i = k; j != k && i < N; ++i)
+ std::swap(G[W*j+i], G[W*k+i]);
+ break;
+ }
+ }
+ // keep searching for suitable column for pivot
+ // beware: this will use columns >= K if necessary.
+ for (int j = k + 1; !G[W*k+k] && j < N; ++j) {
+ for (int h = k; h < K; ++h) {
+ if (G[W*h+j]) {
+ // account column swap
+ std::swap(perm[k], perm[j]);
+ for (int i = 0; i < K; ++i)
+ std::swap(G[W*i+k], G[W*i+j]);
+ for (int i = k; h != k && i < N; ++i)
+ std::swap(G[W*h+i], G[W*k+i]);
+ break;
+ }
+ }
+ }
+ assert(G[W*k+k]);
+ // zero out column entries below pivot
+ for (int j = k + 1; j < K; ++j)
+ if (G[W*j+k])
+ for (int i = k; i < N; ++i)
+ G[W*j+i] ^= G[W*k+i];
+ }
+ }
+ void systematic()
+ {
+ for (int k = K-1; k; --k)
+ for (int j = 0; j < k; ++j)
+ if (G[W*j+k])
+ for (int i = k; i < N; ++i)
+ G[W*j+i] ^= G[W*k+i];
+ }
+ void encode()
+ {
+ for (int i = K; i < N; ++i)
+ codeword[i] = codeword[0] & G[i];
+ for (int j = 1; j < K; ++j)
+ for (int i = K; i < N; ++i)
+ codeword[i] ^= codeword[j] & G[W*j+i];
+ }
+ void flip(int j)
+ {
+ for (int i = 0; i < W; ++i)
+ codeword[i] ^= G[W*j+i];
+ }
+ static int metric(const int8_t *hard, const int8_t *soft)
+ {
+ int sum = 0;
+ for (int i = 0; i < W; ++i)
+ sum += (1 - 2 * hard[i]) * soft[i];
+ return sum;
+ }
+public:
+ void operator()(int *rank, uint8_t *hard, const int8_t *soft, const int8_t *genmat)
+ {
+ for (int i = 0; i < N; ++i)
+ perm[i] = i;
+ for (int i = 0; i < N; ++i)
+ softperm[i] = std::abs(std::max<int8_t>(soft[i], -127));
+ std::sort(perm, perm+N, [this](int a, int b){ return softperm[a] > softperm[b]; });
+ for (int j = 0; j < K; ++j)
+ for (int i = 0; i < N; ++i)
+ G[W*j+i] = genmat[N*j+perm[i]];
+ row_echelon();
+ systematic();
+ for (int i = 0; i < N; ++i)
+ softperm[i] = std::max<int8_t>(soft[perm[i]], -127);
+ for (int i = N; i < W; ++i)
+ softperm[i] = 0;
+ for (int i = 0; i < K; ++i)
+ codeword[i] = softperm[i] < 0;
+ encode();
+ for (int i = 0; i < W; ++i)
+ candidate[i] = codeword[i];
+ score[0] = metric(codeword, softperm);
+ int count = 1;
+ int worst = -1;
+ for (int i = 0; i < 2*L; ++i)
+ scoreperm[i] = i;
+ auto update = [this, &count, &worst]() {
+ int met = metric(codeword, softperm);
+ if (met > worst) {
+ score[scoreperm[count]] = met;
+ for (int i = 0; i < W; ++i)
+ candidate[scoreperm[count]*W+i] = codeword[i];
+ if (++count >= 2*L) {
+ std::nth_element(scoreperm, scoreperm+(L-1), scoreperm+2*L, [this](int a, int b){ return score[a] > score[b]; });
+ worst = score[scoreperm[L-1]];
+ count = L;
+ }
+ }
+ };
+ for (int a = 0; O >= 1 && a < K; ++a) {
+ flip(a);
+ update();
+ for (int b = a + 1; O >= 2 && b < K; ++b) {
+ flip(b);
+ update();
+ for (int c = b + 1; O >= 3 && c < K; ++c) {
+ flip(c);
+ update();
+ for (int d = c + 1; O >= 4 && d < K; ++d) {
+ flip(d);
+ update();
+ for (int e = d + 1; O >= 5 && e < K; ++e) {
+ flip(e);
+ update();
+ for (int f = e + 1; O >= 6 && f < K; ++f) {
+ flip(f);
+ update();
+ flip(f);
+ }
+ flip(e);
+ }
+ flip(d);
+ }
+ flip(c);
+ }
+ flip(b);
+ }
+ flip(a);
+ }
+ std::sort(scoreperm, scoreperm+count, [this](int a, int b){ return score[a] > score[b]; });
+ for (int j = 0, r = 0; j < L; ++j) {
+ if (j > 0 && score[scoreperm[j-1]] != score[scoreperm[j]])
+ ++r;
+ rank[j] = r;
+ for (int i = 0; i < N; ++i)
+ set_be_bit(hard+j*((N+7)/8), perm[i], candidate[scoreperm[j]*W+i]);
+ }
+ }
+};
+
}
diff --git a/tests/osld_regression_test.cc b/tests/osld_regression_test.cc
new file mode 100644
index 0000000..aed2cd5
--- /dev/null
+++ b/tests/osld_regression_test.cc
@@ -0,0 +1,190 @@
+/*
+Regression Test for ordered statistics list decoding
+
+Copyright 2024 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <random>
+#include <cassert>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include "osd.hh"
+#include "galois_field.hh"
+#include "bose_chaudhuri_hocquenghem_encoder.hh"
+#include "bose_chaudhuri_hocquenghem_decoder.hh"
+
+int main()
+{
+ const int L = 16;
+ const bool oracle = true;
+#if 1
+ // BCH(127, 64) T=10
+ const int O = 2;
+ const int N = 127;
+ const int K = 64;
+ const int NR = 20;
+ const int loops = 1000;
+ const double low_SNR = -5;
+ const double high_SNR = 5;
+ const double QEF_SNR = 1.5;
+ typedef CODE::GaloisField<7, 0b10001001, uint8_t> GF;
+ std::initializer_list<int> minpols {
+ 0b10001001, 0b10001111, 0b10011101,
+ 0b11110111, 0b10111111, 0b11010101,
+ 0b10000011, 0b11101111, 0b11001011,
+ };
+#endif
+#if 0
+ // NASA INTRO BCH(15, 5) T=3
+ const int O = 1;
+ const int N = 15;
+ const int K = 5;
+ const int NR = 6;
+ const int loops = 100000;
+ const double low_SNR = -7;
+ const double high_SNR = 7;
+ const double QEF_SNR = 3.5;
+ typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
+ std::initializer_list<int> minpols {
+ 0b10011, 0b11111, 0b00111
+ };
+#endif
+ GF instance;
+ const int NW = (N+7)/8;
+ const int KW = (K+7)/8;
+ const int PW = (N-K+7)/8;
+ int8_t genmat[N*K];
+ CODE::BoseChaudhuriHocquenghemGenerator<N, K>::matrix(genmat, true, minpols);
+ CODE::BoseChaudhuriHocquenghemEncoder<N, K> bchenc(minpols);
+ CODE::LinearEncoder<N, K> linenc;
+ CODE::OrderedStatisticsListDecoder<N, K, O, L> osddec;
+ CODE::BoseChaudhuriHocquenghemDecoder<NR, 1, K, GF> bchdec;
+ GF::value_type erasures[NR];
+ uint8_t message[KW], parity[PW], decoded[NW*L], codeword[NW];
+ int8_t orig[N], noisy[N];
+
+ std::random_device rd;
+ typedef std::default_random_engine generator;
+ typedef std::uniform_int_distribution<int> distribution;
+ auto data = std::bind(distribution(0, 255), generator(rd()));
+
+ int rank[L];
+ double symb[N];
+ double min_SNR = high_SNR;
+ for (double SNR = low_SNR; SNR <= high_SNR; SNR += 0.1) {
+ //double mean_signal = 0;
+ double sigma_signal = 1;
+ double mean_noise = 0;
+ double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
+
+ typedef std::normal_distribution<double> normal;
+ auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
+
+ int awgn_errors = 0;
+ int quantization_erasures = 0;
+ int uncorrected_errors = 0;
+ int ambiguity_erasures = 0;
+ int bchdec_errors = 0;
+ for (int l = 0; l < loops; ++l) {
+ for (int i = 0; i < KW; ++i)
+ message[i] = data();
+
+ linenc(codeword, message, genmat);
+
+ for (int i = 0; i < K; ++i)
+ assert(CODE::get_be_bit(codeword, i) == CODE::get_be_bit(message, i));
+
+ bchenc(message, parity);
+
+ for (int i = K; i < N; ++i)
+ assert(CODE::get_be_bit(codeword, i) == CODE::get_be_bit(parity, i-K));
+
+ for (int i = 0; i < N; ++i)
+ orig[i] = 1 - 2 * CODE::get_be_bit(codeword, i);
+
+ for (int i = 0; i < N; ++i)
+ symb[i] = orig[i];
+
+ for (int i = 0; i < N; ++i)
+ symb[i] += awgn();
+
+ // $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
+ // $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
+ double DIST = 2; // BPSK
+ double fact = DIST / (sigma_noise * sigma_noise);
+ for (int i = 0; i < N; ++i)
+ noisy[i] = std::min<double>(std::max<double>(std::nearbyint(fact * symb[i]), -128), 127);
+
+ osddec(rank, decoded, noisy, genmat);
+ bool unique = rank[0] != rank[1];
+
+ if (oracle) {
+ for (int j = 0; j < L; ++j) {
+ bool found = true;
+ for (int i = 0; i < N; ++i)
+ found &= CODE::get_be_bit(decoded+NW*j, i) == CODE::get_be_bit(codeword, i);
+ if (found) {
+ for (int i = 0; j && i < NW; ++i)
+ decoded[i] = decoded[NW*j+i];
+ unique = true;
+ break;
+ }
+ }
+ }
+
+ for (int i = 0; i < N; ++i)
+ awgn_errors += noisy[i] * orig[i] < 0;
+ for (int i = 0; i < N; ++i)
+ quantization_erasures += !noisy[i];
+ for (int i = 0; i < N; ++i)
+ uncorrected_errors += CODE::get_be_bit(decoded, i) != CODE::get_be_bit(codeword, i);
+ if (!unique)
+ ambiguity_erasures += N;
+
+ for (int i = 0; i < K; ++i)
+ CODE::set_be_bit(message, i, noisy[i] < 0);
+ for (int i = K; i < N; ++i)
+ CODE::set_be_bit(parity, i-K, noisy[i] < 0);
+ int erasures_count = 0;
+ for (int i = 0; erasures_count < NR && i < N; ++i)
+ if (!noisy[i])
+ erasures[erasures_count++] = i;
+
+ bchdec(message, parity, erasures, erasures_count);
+
+ for (int i = 0; i < K; ++i)
+ bchdec_errors += CODE::get_be_bit(message, i) != CODE::get_be_bit(codeword, i);
+ for (int i = K; i < N; ++i)
+ bchdec_errors += CODE::get_be_bit(parity, i-K) != CODE::get_be_bit(codeword, i);
+ }
+
+ double bit_error_rate = (double)uncorrected_errors / (double)(N * loops);
+ if (!uncorrected_errors && !ambiguity_erasures)
+ min_SNR = std::min(min_SNR, SNR);
+
+ int MOD_BITS = 1; // BPSK
+ double code_rate = (double)K / (double)N;
+ double spectral_efficiency = code_rate * MOD_BITS;
+ double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
+
+ double bchdec_ber = (double)bchdec_errors / (double)(N * loops);
+
+ if (0) {
+ std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
+ std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
+ std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
+ std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
+ std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
+ std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
+ std::cerr << bit_error_rate << " bit error rate." << std::endl;
+ std::cerr << bchdec_ber << " BCH decoder bit error rate." << std::endl;
+ } else {
+ std::cout << SNR << " " << bit_error_rate << " " << bchdec_ber << " " << EbN0 << std::endl;
+ }
+ }
+ std::cerr << "QEF at: " << min_SNR << " SNR" << std::endl;
+ assert(min_SNR < QEF_SNR);
+ std::cerr << "Ordered statistics list decoding regression test passed!" << std::endl;
+ return 0;
+}
Served by rngit 1.4.1 - Generated in 0.04s